A new integral representation for quasiperiodic fields and its 
application to two-dimensional band structure calculations 



Alex Barnett 3 *, Leslie Greengard b 

"Department of Mathematics, Dartmouth College, Hanover, NH, 03755, USA 
h Courant Institute, New York University, 251 Mercer St, NY, 10012, USA 



Abstract 

In this paper, we consider band-structure calculations governed by the Helmholtz or Maxwell 
equations in piecewise homogeneous periodic materials. Methods based on boundary integral 
equations are natural in this context, since they discretize the interface alone and can achieve 
high order accuracy in complicated geometries. In order to handle the quasi-periodic conditions 
which are imposed on the unit cell, the free-space Green's function is typically replaced by its 
quasi-periodic cousin. Unfortunately, the quasi-periodic Green's function diverges for families 
of parameter values that correspond to resonances of the empty unit cell. Here, we bypass this 
problem by means of a new integral representation that relies on the free-space Green's function 
alone, adding auxiliary layer potentials on the boundary of the unit cell itself. An important as- 
pect of our method is that by carefully including a few neighboring images, the densities may be 
kept smooth and convergence rapid. This framework results in an integral equation of the second 
kind, avoids spurious resonances, and achieves spectral accuracy. Because of our image struc- 
ture, inclusions which intersect the unit cell walls may be handled easily and automatically. Our 
approach is compatible with fast-multipole acceleration, generalizes easily to three dimensions, 
and avoids the complication of divergent lattice sums. 
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1. Introduction 

A number of problems in wave propagation require the calculation of quasi-periodic solu- 
tions to the governing partial differential equation in the frequency domain. For concreteness, 
let us consider the two-dimensional (locally isotropic) Maxwell equations in what is called TM- 
polarization 1271 12811 . In this case, the Maxwell equations reduce to a scalar Helmholtz equation 



Au(x, y) + cj ep u(x, y) = 0, 



(1) 
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Figure 1 : a) Problem geometry: an infinite dielectric crystal, in the case where the inclusion Q. lies within a parallelogram 
unit cell U. The (shaded) set of all inclusions in the lattice, denoted by Q.\ in the text, has refractive index n, while the 
white region has index 1. b) Sketch of our quasi-periodizing scheme: we make use of layer potentials on the left (L) and 
bottom (B) walls, extended to the additional segments shown, which form a skewed 'tic-tac-toe' board, as well as the 
near neighbor images of fl, outlined in solid lines. 



where e and fi are the permittivity and permeability of the medium, respectively, and we have 
assumed a time dependence of e~ ,aJt at frequency to > 0. Given a solution u to (Q}, it is straight- 
forward to verify that the corresponding electric and magnetic fields E, H of the form 

E(x,y,z) = E(:c,y) = (0,0, u(x,y)) 

H(x,y,z) = H(x,y) = - — (u y (x,y) ) -u x (x,y),0) 



satisfy the full system 



VxE 
VxH 



—icoeE . 



We are particularly concerned with doubly periodic materials whose refractive index n = yfql is 
piecewise constant (Fig. fl}. Such structures are typical in solid state physics, and are of particular 
interest at present because of the potential utility of photonic crystals, where the obstacles are 
dielectric inclusions with a periodicity on the scale of the wavelength of light [28]. Photonic 
crystals allow for the control of optical wave propagation in ways impossible in homogeneous 
media, and are finding a growing range of exciting applications to optical devices, filters [21], 
sensors, negative-index and meta-materials jj^l . and solar cells 10]. 

We assume that the crystal consists of a periodic array of obstacles (£2a) with refractive index 
n + 1, embedded in a background material with refractive index n = 1 (denoted by R 2 \ Qa). We 
then rewrite (fl} as a system of Helmholtz equations 



(A + n 2 co 2 )u = 



in Qa 



(A + ar)u 







in . 



(2) 
(3) 



The expression Qa, above, is used to denote the closure of the domain Qa (the union of the 
domain and its boundary 3Qa). In this formulation, we must also specify conditions at the 
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material interfaces. These are derived from the req uired continuity of the tangential components 
of the electric and magnetic fields across <9Qa B27L 12811 . yielding 

u, u„ continuous across 8Q,\ (4) 

where u„ = du/dn is the outward-pointing normal derivative. 

The essential feature of doubly periodic microstructures in 2D (or triply periodic microstruc- 
tures in 3D) is that, at each frequency, there may exist traveling wave solutions (Bloch waves) 
propagating in some direction defined by a vector k. 

Definition 1. Bloch waves are nontrivial solutions to ©-(HI that are quasiperiodic, in the sense 
that 

u(x) = e**u(x) , (5) 

where li is periodic with the lattice period andk — (k x , k y ) is real-valued, k is referred to as the 
Bloch wavevector. 



Bloch waves characterize the bulk optical properties at frequency oj; they are analogous to 
plane waves for free space. If such waves are absent for all directions k for a given a>, then the 
material is said to have a band-gap 14811 ). The size of a band-gap is the length of the frequency 
interval [a»i, 0J2] in which Bloch waves are absent. Crystal structures with a large band-gap are 
'optical insulators' in which defects may be used as guides [28], with the potential for enabling 
high-speed integrated optical computing and signal processing. 

Definition 2. The band-structure of a given crystal geometry is the set of parameter pairs (u>,k) 
for which nontrivial Bloch waves exist. 

The numerical prediction of band structure is a computationally challenging task, yet essen- 
tial to the design and optimization of practical devices. It requires characterizing the nontrivial 
solutions to a homogeneous system of partial differential equations (0, ([3]) subject to homoge- 
neous interface and periodicity conditions ©, ((5]) in complicated geometry. Solving this eigen- 
value problem is the focus of our paper. 

In the next section, we briefly review existing approaches, and in section [3] we present and 
test a method that relies on the quasi-periodic Green's function. We introduce our new mathe- 
matical formulation in section [4] Numerical results are presented in section [3] and we conclude 
in section [6] with some remarks about the potential for wider application of this approach. 



2. Existing approaches 

In order to pose the band-structure problem as an eigenvalue problem on the unit cell U (see 
Fig. Q]), we will require some additional notation. The nonparallel vectors ei,e2 e IR 2 define 
a Bravais lattice A := {mei + ne.2 : m,n e Z}. Given a smooth, simply connected inclusion 
Q c R 2 , we may formally define the corresponding dielectric crystal by Qa := {Q. + d : d e A). 
As indicated above, we assume that Qa has refractive index n + 1, and that the background 
M 2 \ Qa has refractive index 1. For the moment, we assume that Q c U as illustrated in Fig.Q] 
We will discuss the case of Q crossing dU in Section I5TI 
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The quasi-periodicity condition © can be rewritten as a set of boundary conditions on the 
unit cell U, coupling the solution on the left (L) and right (L + ei) walls, as well as on the bottom 
(B) and top (B + £2) walls. More precisely, if we define 



a:=k-ei, a := e' a , b := k • e 2 , B := e' b , 
then quasi-periodicity is written 

M| L+e , = au\ L (6) 

"nli+e, = OlU n \ L (7) 

«lB+e 2 = /Mb (8) 

U n \B+e 2 = Bu n \ B , (9) 

where the normals have the senses shown in Fig.Q] 

The homogeneous equations (f2]i-(|4]i, (If)}-© define a partial differential equation (PDE) eigen- 
value problem on the torus U. By convention, the band structure or Bloch eigenvalues are gener- 
ally defined as the subset of the parameter space {(«, a,b) : cd > 0,—n < a < n,-n < b < n] for 
which nontrivial solutions u : U — > C exist. The earlier definition of band-structure, based on 
(0, allows for arbitrary values of k. It is clear, however, that one only needs to consider a single 
period of k's projection onto ei , e2, which we have denoted by a, b, to characterize the entire set 
(w, k) of nontrivial Bloch waves. This domain {(a,b) : —n < a < n, -n < b < jt} is (essentially) 
what is referred to as the Brillouin zone. 

Because the PDE is elliptic and U is compact, for each k there is a discrete set of eigenvalues 
{<y ; (k)}°lj, counting multiplicity, accumulating only at infinity. Each a»/k) is continuous in k, 
so that the bands form sheets. 

Popular numerical methods for band structure calculations are reviewed in ll28"ll . Broadly 
speaking, they may be classified as either time-domain or frequency domain schemes. In the first 
case, an initial pulse is evolved via the full wave equation (typically using a finite-difference or 
finite-element approximation). If the simulation is sufficiently long, Fourier transformation in 
the time variable then reveals the full band structure. In the second case, the eigenvalue problem 
©-([Hi, (HJ-© is discretized directly. Such frequency domain schemes can be further categorized 
as 

1 . PDE-based methods, which involve discretizing the unit cell using finite difference or finite 
element methods Ji, 19, 20ll . 



2. plane-wave methods which expand the function u in (0 as a Fourier series, and apply the 
partial differential operator in Fourier space 128112911 . 

3. semi-analytic multipole expansion methods which apply largely to cylindrical or spherical 
inclusions fiol 431. 



methods which use a basis of particular solutions to the PDE at a given frequency a> and 
enforce both interface and boundary conditions as a linear system, such as the "multiple 



multipole" or "transfer-matrix" method II23L 14611 . and 



5. boundary integral (boundary element) methods 04911 . which includes the method described 
here. 

For a fixed k, methods of type (1) and (2) result in large, sparse generalized eigenvalue prob- 
lems whose lowest few eigenvalues approximate the first few bands a>j(k). They have the advan- 
tage that they couple easily to existing robust linear algebraic techniques. PDE-based methods, 
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however, require discretization of the entire cell in a manner that accurately resolves the geom- 
etry of the inclusion £1 Plane-wave methods, which perform extremely well when the index 
of refraction n is smooth, have low order convergence when n is piecewise constant, as in the 
present setting. Both require a large number of degrees of freedom. 

Methods of type (3), (4) or (5), on the other hand, represent the solution using specialized 
functions (solutions of the PDE) whose dependence on at is nonlinear. As a result, they can 
be much more efficient and high-order accurate, dramatically reducing the number of degrees 
of freedom required. Unfortunately, however, they result in a nonlinear eigenvalue problem 
involving all the parameters a>, a and b, and somewhat non-standard techniques are required to 
find values of the parameters for which the system of equations is singular M47I1 . 

We are particularly interested in using boundary integral methods (BIEs), since they eas- 
ily handle jumps in the index in complicated geometry, have a well understood mathematical 
foundation, and can achieve rapid convergence, limited only by the order of accuracy of the 
quadrature rules used. High order accuracy is important, not only because of the reduction in the 
size of the discretized problem, but in carrying out subsequent tasks, such as sensitivity analyses 
ifTrh through the numerical approximation of derivatives, and the computation of band slopes 
(group velocity), and band curvatures (group dispersion). 

There is surprisingly little historical literature on using BIE for band structure calculations, 
although the last few years have begun to see some activity in this direction (see, for example, 
Ugh). There is, however, an extensive literature on integral equations for scattering from periodic 
structures, which we do not seek to review here. For some recent work and additional references, 
see QUI. 



3. Integral equations based on the quasi-periodic Green's function 

An elegant approach to designing integral representations for quasiperiodic fields involves 
the construction of the Green's function that imposes the desired conditions ©-(O exactly. We 
first need some definitions lfl6l 41]. At wavenumber a> > 0, the free space Green's function for 
the Helmholtz equation, G is defined by -(A + <d 2 )G = do where So is the Dirac delta function 
centered at the origin. In 2D, this yields 



G(x) = G M (x) = -h£\oj\x\), x e M 2 \ {0}, (10) 



where H^ 1 is the outgoing Hankel function of order zero. By formally summing over images of 
the Green's function placed on the lattice A, with correctly assigned phases, we get an explicit 
expression for the quasi-periodic Greens function 

G QP (x) = Yj e ik d G(x - d) = Yj a'TG(x - me, - ne 2 ) . (11) 

dsA m,neZ 

We leave it to the reader to verify that G QP does, indeed, satisfy (O-®. One small caveat: the 
series in (fTTb is conditionally convergent for real u). The physically meaningful limit is taken 
by assuming some dissipation a> = co + is in the limit e — > + (see 11811 for a more detailed 
discussion). It will be useful to distinguish between the copy of the Green's function sitting 
in the unit cell U and the set of all other images. For this, we define the "regular" part of the 
quasi-periodic Green's function by 

G£ p (x) = Y a'T^x-md -ne 2 ). (12) 



This function is a smooth solution to the Helmholtz equation within U and clearly satisfies 



G QP (x) = G(x) + G r Qe (x) . (13) 

A spectral representation also exists J^,[T§], built from the plane-wave eigenfunctions of the 
quasi-periodic torus U : 

I gi(k+q)x 

Vol(t/) ^ |k + q| 2 - w 2 

Here, A* := \mr\ +nr2 : m, « e Z] is the reciprocal lattice with vectors r,- defined by e,-ry = Indij 
for z, 7 = 1,2. From the denominators in (fT4b it is clear that G QP may blow up for specific 
combinations of a> and k. The quasiperiodic Green's function is, in fact, well-defined if and only 
if those parameters satisfy the following non-resonance condition. 

Definition 3 (empty resonance). A parameter set (cj,k), equivalently {u,a,b), is empty reso- 
nant if (±> = \k + q\for some q € A*, otherwise it is empty non-resonant. 

Our terminology comes from the fact that the blow-up in G QP is physically the resonance of the 
'empty' unit cell U, with refractive index 1 everywhere and quasi-periodic boundary conditions. 
That is, G QP is undefined if and only if (w, a, b) lies on the band structure of the empty unit 
cell. The blow-up of the Green's function is less apparent from (fTTT >. but is manifested in the 
divergence of the series, even in the limit co = a> + is with e — > + . 

It will be convenient sometimes to refer to a Green's function as a function of two variables, 
with G(x, y) := G(x-y), and G QP (x,y) := G QP (x-y). Then, for each y e R 2 , the function G QP (-,y) 
is quasi-periodic. 

We now represent solutions to the PDE eigenvalue problem (f2|i-(|4|i, ©-fig} by the layer po- 
tentials, 

S (n0J) cr + D (,m) T in a 
S QP 'cr + D qp 't in U \ D. 

where the usual single and double layer densities [16] at any wavenumber oj > are defined by 

(5 (w V)(x) = f G^(x,y)o-(y)ds y (16) 
Jan 



r dG (u>> 

__(x,y)T(y)^ y (17) 
Jan oriy 



and their quasi-periodized versions are likewise 

(S$cr)(x) = f G^(x,yMy)^ y (18) 
Jan 

(D ( ^T)(X) = -^( X) y) T ( y )^ (19) 

Jan dn y 

Here ds is the usual arc length measure on <9Q, and the derivatives are with respect to the second 
variable in the outward surface normal direction at y. It is clear [ 16] that the above four fields 
satisfy the Helmholtz equation at wavenumber to in both Q and U \ Q. Note that we have 
chosen a non-periodized representation within the inclusion O in ( fT31 ). which has some analytic 
advantages (see Theorem[4]and the last paragraph in the Appendix). 
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Since u in dT5l > satisfies (O, ([3J, and all that remains is to solve for densities cr, r such 

that the matching conditions (0]l are satisfied, which we now address. 

Using superscripts + and - to denote limiting values on dQ, approaching from the positive 
and negative normal side respectively, we use the field ( fTST l and the standard jump relations for 
single and double layer potentials [ 16|,|22|] to write 



u + - u 




/ 







u n — u n 


■( 





/ 


+ 



D ^)_ D (noj) S<- naj) -S Qf 



QP 
(to) ' 



1 QP 1 1 ^QP 



) 


T 




-O" 



=: A QP ?7 



(20) 



Here / is the identity operator, while S and D are defined to be the limiting boundary integral 
operators (maps from C(dCl) — > C(<90)) with the kernels S and £) interpreted in the principal 
value sense. (S is actually weakly singular so the limit is already well defined. A standard 
calculation | Tg, 2^1 shows that D is weakly singular as well). The hypersingular operator T has 
the kernel gn gn y and is unbounded as a map from C(<9Q) — > C{dQ). In these definitions, as in 
(IT6l-(fT9]l. it is implied that G inherits the appropriate superscripts and subscripts from S , D and 
T. Finally, * indicates the adjoint. The amounts by which the material matching conditions fail 
to be satisfied, 



it — it 

ll n Ut, 



(21) 



is a column vector of functions which we call the mismatch. We summarize the linear system 
(f20t by m = A QP j] where rj :- [r; -cr]. It is important to note that the difference of hypersingular 
kernels, T$ - T (n "\ in ® is only weakly singular lfl6l Sec. 3.8]. This cancellation, achieved 
here by using the same pair of densities inside as outside the inclusion, is well known ll44ll . 
The result is that A QP is a compact perturbation of the identity and (l20b is a Fredholm system of 
integral equations of the second kind. 

In the above scheme, we might hope that if it is possible to find nontrivial densities r\ whose 
field u gives zero mismatch m for a set of parameters (at, a, b), then that set is a Bloch eigenvalue. 
Indeed (as with the case of simpler domain eigenvalue problems lf39L Sec. 8]) we have a stronger 
result. 



Theorem 4. Let (10, a,b) be empty non-resonant. Then (a>, a, b) is a Bloch eigenvalue if and only 
?/NullA QP + {0} . 



The proof occupies Appendix | Appendix A| This suggests the core of a numerical scheme: at 
each of a sampling (e.g. a grid) of parameters (w, a, b), find the lowest singular value cr min (A QP ) of 
a matrix discretization A QP of A QP . The band structure will then be found where <r min (A QP ) is close 
to zero. 



3.1. Discretization of the integral operators 

Since the goal of this work is to explore periodization, we limit ourselves to the simplest case 
of dQ being smooth. The methods of this paper extend without much effort to other shapes, but 
the quadrature issues become more involved. Recalling ( fT3l , note that the kernels in (l20b are the 
sum of a component due to G which is weakly singular, plus the remainder due to G' Qf which is 
smooth (analytic). We will make use of a Nystrom discretization using the spectral quadrature 
scheme of Kress lf3lll for G and the trapezoidal rule for G£, p . 
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We first remind the reader of the periodic trapezoidal Nystrom scheme lf33ll . in the context of 
a general second kind boundary integral equation 



yu(x) + I k(x, y)(x(y)ds y = f(x), x e <9Q, 
Jdn 



where d£l is parametrized by the 27r-periodic analytic function z : [0, 2n) — > R 2 . Changing 
variable gives 

pis) + I K(s, t)n(t)dt = f(s), s G [0, 2tt), 



I 

Jo 

where K(s,t) :- k(z(s),z(t))\z'(t)\ and z' = dz/dt. Choosing N quadrature points tj = 2nj/N 
with equal weights 2n/N gives the N-by-N linear system for the unknowns //. , which approxi- 
mate the exact values p{tj), as 

? N 

+ f X K(tk > '^T = f (tk) > k=\,...,N. (22) 



By Anselone's theory of collectively compact operators 1 33 1 , the convergence of errors \f/p - 

fJ.(tj)\ inherits the order of the quadrature scheme applied to the exact integrand K(s, which is 
analytic when k and / are. 

Remark 5. For analytic integrands, the periodic trapezoidal rule has exponential convergence 
with error 0(e~ 2yN ) where y is the smallest distance from the real axis of any singularity in the 
analytic continuation of the integrand. [33, Thm. 12.6]. 

The above discretization is used to populate the matrix entries in (f20b that are due to the 
smooth compoment G r QP . (We explain how to compute this kernel itself in Section [3~2l ) 

For non-smooth kernels, such as G, the rule ( l22b must be replaced by a quadrature that cor- 
rectly accounts for the singularity in order to retain high order accuracy. There are a variety 
of such schemes, such as those of J2, HH [3(1 • By fixing the order of accuracy, they allow for 
straightforward coupling to fast multipole acceleration Jl2 . 13. 14. 42j] by making local mod- 
ifications of a simple underlying quadrature rule (such as the trapezoidal rule or a composite 
Gaussian rule). In the present context, we ignore such considerations and use a global rule due 
to Kress 113 111 that achieves spectral accuracy in the logarithmically singular case. 

The essential idea of Kress' scheme (after transformation of variables to the interval [0, 2n]) 
is to split a logarithmically singular kernel K(s, t) in the form 

K(s,t) = log ^4 sin 2 ^)^i(i,f) + ^2(s,0 (23) 

with K\ and K2 periodic and analytic. K2 is (again) handled with the trapezoidal rule. For K\, 
the Kussmaul-Martensen quadrature rule is spectrally accurate: 

JT log(4sin 2 ^)g(f)A *Y^Rf\s)g(tj) (24) 

with quadrature weights (deriving from the Fourier series of the log factor) given by 

W2-1 2 2 N 

R ( P(s) = - V -cosm^-O) - - cos -(s-tj). (25) 
J ^ m N 2 



a) b) c) 




Figure 2: Convergence for quasi-periodic Greens function scheme of Sec. [3] a) Absolute error in o" min (AQ P ) vs N the num- 
ber of quadrature nodes on dil, for Bloch parameters a = n/2 and the two different b values labeled. The unit cell with 
ei = (1, 0), e2 = (0.4, 1), and inclusion, described by the radial function r(8) = 0.2(1 + 0.3 cos 38), are shown in the inset. 
Index is n = 3 and frequency a/ = 4.5. For b = 2, error is taken relative to the converged value 0.01879908530381247; 
for b = bo, relative to 0. The matrix Aq P has 2-norm of about 25. b) o" min (AQ P ) vs difference in parameter b from the 
Bloch eigenvalue bg, for several different numbers of quadrature points N. Note the horizontal log scale. This shows that 
it is the convergence rate at the Bloch eigenvalue that controls the accuracy with which the minimum can be found, c) 
Relative error (+ symbols) in evaluation of lattice sum S3 by the method of Sec. l3.2l vs the maximum order L in )27t . 
Parameters are as in Table II of HI, whose claim S3 = 2.13097899279352 + 5.66537068305984; is taken as the true 
value. Also, relative error (o symbols) for § 3 which excludes the 3x3 block of neighbors (parameters are the same; true 
value is taken as the converged value at L = 50). 



Thus, the matrix elements in discretizing (l23l are Kit^, tj) = Ry_ k ^0)Ki(tk, tj)+K2(tk, tj). Finally, 
it is always the difference of two hypersingular operators T that appears in the integral equation 
(120) . This difference is only logarithmically singular, so that Kress' rule can be used for every 
block of (EDI . We refer the reader to 113 ill for further details. 

In summary, a matrix discretization A QP of A QP is formed by using the above quadrature rules 
for each of the 2-by-2 integral operator blocks in (1201 1. This matrix maps density values to field 
values. However, in order to create a matrix whose singular values approximate those of A QP 
we must instead normalize such that 2A^-dimensional Euclidean 2-norms correctly approximate 
L 2 (dQ)-norms. This is done by symmetrizing using quadrature weights to give our final matrix 

A QP = W l/2 A QP W- l/2 (26) 

where W is diagonal with diagonal elements Wj = wj+n = (2n/N)\z'(tj)\, for j = l,...,N. 

The net result of the preceding discussion is that with the use of specialized quadratures on 
smooth boundaries, the singular values of A QP are spectrally accurate approximations to those 
of A QP . We demonstrate this convergence for a small trefoil-shaped inclusion in Fig. [2^; the 
convergence is spectral, until the error is approximately machine precision times the matrix 2- 
norm. The rate appears to be faster at a Bloch eigenvalue (in this case on the fourth band) than 
far from one. Fig. EJ3 shows that the minimum locates the parameter b to 14 digit accuracy for 
JV>70. ^ 
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3.2. New method for evaluation of the quasi-periodic Greens function 

In order to compute the elements of A QP , one must evaluate G' QP defined by (fT2b ; in this 
section, we present a surprisingly simple (and apparently new) method for this. Since the sums 
(fTTT > and (TBI converge too slowly to be numerically useful, many sophisticated schemes have 
been devised. Some of these are based on the Fourier representation (such as [9]), but most are 
based on the observation that 



G' Qf (r,6)=Y j S l J l (ajry w . 



(27) 



l=-L 



where (r, 0) are the usual polar coordinates, and 7/ the regular Bessel function of order /. As 
L — > oo, this expression is uniformly convergent in the unit cell U, as long as there exists a circle 
about the origin which contains U but encloses no points in A \ {0}. The coefficients 5/ in this 
expansion are know as lattice sums, given by 

S l = J] a m r3"Hf\u>r mn )e- iW >"", 

(m,«)#(0,0) 



where (r mn , 6 mn ) are the polar coordinates of mei + nea, and Hj' is the outgoing Hankel function 
of order /. Thus, the issue of evaluating G r Ql , has been reduced to that of tabulating the lattice 
sums. This problem itself has a substantial literature (see, for example, 1 15, 18 , 34, 38, 40ll ). 
Nevertheless, very few papers discuss the problem of empty resonances, at which point the lat- 
tice sums S i blow up. One notable exception is the work of Linton and Thompson 1I35I1 . who 
analyze this blowup for periodic one-dimensional arrays in two dimensional scattering. They 
also propose a regularization method to overcome it. 

We present here the construction of a small linear system whose solution yields the lattice 
sums rather easily (away from empty resonances). In physical terms, we compute the field 
induced by the free-space Green's function G, determine how it fails to satisfy quasi-periodicity, 
and use the representation ( IZTb to enforce quasi-periodicity numerically. More precisely, given a 
field m, we define the discrepancy by 



a l u\ L+ei 

- <Z -1 «nlz,+e, 



u\ L - 
u„\l 

u\b - P 1 w|fi+e 2 



U n \B 



U n \B+e 2 



(28) 



We can interpret /, /' as functions on wall L and g, g' as functions on wall B. We construct a 
4M-component column vector d by sampling these four functions at Gaussian quadrature points 



{y { m)Z=\ on L ' and itm ; )„=i on B. If we let the field u(x) = G(x), then for m = 1, . . . , M, the mth 



JBhM 



element of d is G(y^) - aT 1 G(y^ - ' + ei). The remaining 3M entries in d are computed in the 
analogous fashion. 

Now let H be a (complex) matrix of size AM x (2L + 1 ), defined as follows. For / = — L, . . . , L, 
fill the (I + L + l)th column in the same manner as d, but using the field m(x) = Ji(ajr)e' w . Letting 
s := {Si\f_ _ L , it is straightforward to verify that the linear system 



Hs = -d 



(29) 
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yields values for the lattice sums that annihilate the discrepancy induced by the source G. We 
solve the linear system in the least squares sense. This has to be done with some care, since the 
Bessel functions J\ become exponentially small for large I. A simple fix is to right-precondition 
the system by scaling the (/ + L + l)th column of H by the factor p; := l/Ji(mm[a>R, /]), where 
R :- max xe ;y |x| is the unit cell radius. The entire procedure may be interpreted as finding the 
representation ( l27T i which minimizes the L -norm of the discrepancy of the resulting Gqp. 

Fig. [2J3 shows that the error in evaluating S /, for / = 3, has exponential convergence in L. We 
fixed M — 24 (large enough that further increase had no effect). 14 digits of relative accuracy 



are achieved for L > 46, comparable in accuracy to 113 811 . Although the maximum achievable 
accuracy for Si deteriorates exponentially as |/| increases, the resulting accuracy of G' QP computed 
via (l27l > is close to 14 digits everywhere in t/Q We do not claim that our method is optimal in 
terms of speed (although at 0.05 sec to solve for all S / values, it is adequate), merely that it is 
accurate, convenient and robust. To our knowledge it has not been proposed in the literature. 

The convergence rate in the boundary L 2 -norm of expansions such as (l27T i depends on the 
(conformal) distance from the domain to the nearest field singularity (a result of Vekua's theory 
and approximation in the complex plane [8, Ch. 6]). Thus, the rate may be improved by increas- 
ing this distance by removing the rest of the 3x3 block of nearest neighbors from the lattice 
sum, and representing 

L 

G r QP (x) := Yj a^Gix - M - ke 2 ) = J] §iJi(cor)e m . (30) 

O',jt)eZ 2 \(-l,0,l) 2 l=-L 

To solve for {§ /}, the right-hand side of the linear system is now chosen to be the direct summa- 
tion of these neighbors, m(x) = G(x) := 2i*e{-i,Q,i) a'/3 k G(x - jei - k&i)- We may then evaluate 
G QP = G + Gg P . As Fig.[2j: shows, the convergence rate for Si, and hence for G QP , is now a factor 
2-3 better. Hence we use this method below, fixing L = 30. 



3.3. The empty resonance problem 

Given a photonic crystal (inclusion Q. with index «), using the methods of Sections [3~Tl and [3T2l 
we are able to construct the matrix A QP for any given frequency and Bloch parameters {u, a, b). 
Fig. [3^ shows the minumum singular value of this matrix as a function over the (b, of) plane, for 
constant a: the band structure is visible as the zeros of this function. We have also superimposed 
the band structure of the empty unit cell (dotted lines). Theorem|4]guarantees that, away from the 
empty unit cell band structure, no spurious modes will be found, and that no modes are missed. 

However, zooming in to one of the many intersections of the two sets of curves (Fig. [3j)), 
we see that in the neighborhood of the empty band structure, the desired singular values take 
on arbitrary fluctuating values that obscure the theoretical behavior near their intersection. This 
prevents any attempt to locate the desired zero set to an accuracy better than 0( yfsZD, where 
e„ md , is the machine precision. As Fig. [3p shows, this is explained by the blowup of the entries 
of the matrix A QP as one approaches the empty band structure. This, in turn, causes unbounded 
roundoff error when computing small singular values in finite-precision arithmetic. 



'This is to be expected from arguments similar to [4| Eq. (5)]: the residual of the linear system, around 10~ 14 , 
approximates the boundary error norm, which in turn controls the interior error norm when using a basis of particular 
solutions to the Helmholtz equation. 
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Figure 3: Breakdown of quasi-periodic Greens function scheme, for the system of Fig.[2ji except with £2 = (0.5, 1). a) 
Minimum singular value of Aq P vs b = k ■ e2 and w, as a log density plot over a slice with fixed a = 0.8. Dark curves 
indicate the band structure, and superimposed dotted lines the 'empty' band structure where Gqp blows up. b) Zooming 
in by a factor of 10 7 to the region shown by the dot in a), showing failure to resolve band structure at the intersection, 
c) log 1-norm of the matrix Aq P plotted over the same region as b); it is of order the inverse of the distance to the empty 
band structure. 

Remark 6. The above demonstrates a fundamental flaw inherent in the use of the quasi-periodic 
Greens function in band structure problems; there are empty-resonant parameter sets (sheets in 
the space (a>, a, b)) where the desired band structure cannot be computed. Furthermore, loss of 
accuracy is inevitable near these parameter sets. 

This motivates the development of a more robust scheme. 



4. Periodizing using auxiliary densities on the unit cell walls 

4.1. Inclusion images and a new linear system 



Section I3T2I illustrated the fact, well known in the fast multipole literature ll6l ll4ill2lll3Lll8ll . 
that summing the nearest neighbors directly (i.e. excluding them from the quasi-periodic field 
representation) results in much improved convergence rates. This motivates defining general- 
izations of ( fTol l and (TPTT i that include summation over the appropriately phased 3x3 nearest 
neighbor images, as shown in Fig.QJ, 



(<S ( "V)(x) = f J] ffW w) (x,y + jei+fce 2 )o-(y)^ y (31) 
^ an j,te{-l,0,l) 

y a jpk (x ^ + M+kt2)T(y)d (32) 



We now choose a layer potential representation for u that involves only free space kernels: 

J S ( ' m) o- + £> (nu) T inO 

The auxiliary field u QP will be represented by a new set of layer potentials that lie on the "tic-tac- 
toe" stencil of Fig.QJ), consisting of the boundary of U and its closest extensions, none lying in 
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the interior of f/. We will return to this in section Ek2l For the moment, let us denote the unknown 
densities that determine u QP by By construction, the representation d33l satisfies (f2]i and (01 
in U, so that it remains only to impose both the matching/continuity conditions © and quasi- 
periodicity ©-(O. Imposing the mismatch m defined by ( |2TI ) and the discrepancy d defined by 
(1281 on m, the unknowns in (l33l must satisfy a linear system of the form: 

(34) 

where, as before, rj := [r; — cr], We will describe the operators A, B, C, and Q in more detail 
shortly. For the moment, note that if there exists a density [77; which generates a nontrivial 
field with vanishing mismatch and discrepancy, then it is a solution to ©-(O and ©-(O and the 
corresponding parameters (to, a, b) must be a Bloch eigenvalue. Numerical evidence supports the 
following stronger claim, the analog of Theorem[4] 

Conjecture 7. (lo, a, b) is a Bloch eigenvalue if and only j/Null E + {0}. 

This suggests, as in Section [3j computing the band structure by locating the parameter families 
where (a discretization of) E is singular. The point of the new scheme is that it should be robust; 
if the conjecture holds, then (in contrast to the quasiperiodic Green's function approach), there 
will be no spurious parameter values where the method breaks down. 

To discuss the operators in E, we need some additional notation. We assume that the wavenum- 
ber co and quasiperiodicity parameters (a, b) are given. Let W be a curve in M 2 on which single 
and double layer densities are defined, with the corresponding potentials written as 



(S w o-)(x) = f G(x,y)o-(y)ds y (35) 
Jw 

(D w t)(x) = f ^-(x,y)T(y)ds y . (36) 
Jw on y 

Letting V be a (possibly distinct) target curve in R 2 , we define the operators 

(Sv,wo-)(x) = f G(x,y)o-(y)ds y xeV (37) 
Jw 



r dG_ 

Jw dn y 



(ZW)(x) = — (x,y)T(y)J Sy xeV (38) 



y 

y 



(D* vw o-)(x) = f ^(x,yMy)^ y xe V .39s 
Jw Cn x 

(7V,wrXx) = f -^-(x,y)T(y)^ y xeV. (40) 
J w dn x dn y 

When V = W, these operators are to be understood in the principal value sense. By analogy with 
(fJTJ, d32b , versions of these operators whose kernels include the phased sum over 3x3 images 
of the source are indicated with a tilde (~): that is, § v,w, Dy,w> D* vw , and Ty tW - 

We are now in a position to provide explicit expressions for the operators A, B, C, Q in ( l34l . 
Comparing d33l to (15[ , it is clear that the operator A is the same as A QP in < f20b but with the re- 
placement of S%\ and T$\ by S aci.dci, Dan,an and Tan^n, respectively. It is straightforward 
to verify that A is a compact perturbation of the identity. 

13 



Figure 4: Discrepancy cancellation due to neighbor image sums. Each arrow represents the influence of a source density 
on a target segment, a) For the four upper sub-blocks of C, the six nearest source images (dotted) contribute to the 
discrepancy on the left wall L. b) The six nearest source images (dotted) contribute exactly the same field (suitable 
phased) to the right wall (L + ei). The net result is that only the distant sources, shown in bold, have a non-zero effect. 
The same holds for all four upper sub-blocks of C. A rotated version applies to the lower four sub-blocks of C. c,d) 
Contributions to the sub-blocks Qn and gz.s of Q. The seven indicated terms (dotted source segments) cancel in the 
two diagrams, leaving only the contributions from distant wall segments shown in bold. A rotated version applies to the 
sub-blocks Qbl and Q m . 

The operator C describes the effect of the inclusion densities on the discrepancy d. Its eight 
sub-blocks are found by inserting (FJTT l and (l32l i into (1331 ) then evaluating (l28b . giving 



C = 





a 




-S_L,dn+ a 


S L+e u an 
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TL+e u dn 


-&ilKX +a ~ 


1 rj* 
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1 DB+e 2 ,dQ 
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S B+e 2 ,dSl 


. T B ,an- 








l D* 



Consider now the any of the four upper sub-blocks of C. There are nine phased copies of <9fl 
which contribute to the field on the left (L) and right (L+ei) wall. From symmetry and translation 
invariance considerations, however, it is easy to check that the contributions from the six left- 
most images on L (dotted curves in Fig. |4^) are equal to the contributions of the six right-most 
images on L + ei (dotted curves in Fig. [4}}). In the (1, 1) sub-block, for example, we have: 

fa={-l,0,l} 

A rotated version of the analysis applies to the lower four sub-blocks in C. The result is that the 
entries in C involve only source-target interactions at distances greater than the size of the unit 
cell, ensuring the rapid convergence of a representation in terms of smooth functions. 

We next discuss the representation of £ and u QP [^] in more detail, which will determine the 
form of blocks Q and B of the full system matrix E. 

4.2. Choice of auxiliary densities and their images 

The auxiliary field m qp is determined by the choice of layer potentials on the boundary of (and 
outside of) U. We will use double and single layer densities on both the left (L) and bottom (B) 
boundaries of U, as well as on the other segments of the "tic-tac-toe" board in Fig. QJ. More 
precisely, we define the vector of unknowns £ by £ := [tl\ —ctl, t b ; -o~b], and set 

"op = aJ /^{^L+je 1 +kt2 Cr L + £>L+je i+ ke 2 TL) + ^ {^B+je l -t*t 2 0'B + ft+je, +ke 2 T B ) 

fce{-l,0,l) fe(0,l) 

(41) 

The inclusion of the image segments leads to cancellations that are numerically advantageous in 
the operator Q, just as we found that images helped with the operator C in the preceding section. 
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We should first clarify the definition (l28t of the discrepancy functions: field values should 
be interpreted as their limiting values on the wall approaching from inside U, since it is the 
field in U that (l33l l and (PiTt represent. For example, / := u + \i - a~ l u~\i +ei , where, as before 
« ± (x) := lim e ^o+ u(x + en), and n is the normal at x. 

Recall now that the operator Q expresses the effect of the four densities in £ on the four dis- 
crepancy functions /, /', g, g'. If (l4TT i contained only the terms j = k — 0, this would correspond 
to densities cr L and tl placed on L, and cr B and tb placed on B. While this is mathematically 
acceptable, it results in various complicated self-interactions and interactions between segments 
that share a common corner. This would lead to singularities in densities requiring more compli- 
cated discretization and quadrature. Although there has been significant progress in this direction 
(see, for example, 111 ll 12511 '). in the present context we have the luxury of including the ten addi- 
tional image segments in (HTt . which cancel both the self and near-field corner interactions. As 
a result, our implementation is simpler and involves fewer degrees of freedom. The cancellation 
mechanism is shown in Fig. [4] The effect on u + \i of the seven segments touching L, for example, 
cancels the effect on M~|z.+e, of the seven segments touching L + ei, leaving only ten far field 
contributions. 

It is important to note that the local terms due to the jump relations do not cancel: e.g. a 
density function tl placed on L contributes a term 517, to u + \i, while otl placed on L + e\ 
contributes - jcrrz, to M~|z,+ ei ■ These two terms add to contribute tl to /. One may check in this 
fashion that the jump relations contribute an identity to the diagonal sub-blocks of Q. This yields 
the crucial result that Q is the identity plus a compact operator, with the compact part generated 
by interactions at a distance greater than the size of the unit cell. After the above cancellations 
and simplification, we have, 

q _ j + Qll Qlb 
[Qbl Qbb 

where 



Qll = 



2^ ja j /3 k D LMjf . l+kf . 2 

J€{-l,l},foE{-l,0,l) 

^ ja J p k T LMjei+k e 2 
jei-l, L},fee{— 1,0,1} 



- 2u j aJ P kS L,L+M+kei 
j'e{-l,l},Jfce{-l,0,l) 

- j aj P kD *LMM + ke 2 
;£{-l,l},fe(-l,0,l) 



Q 



LB 



/ \ /^(aD^B+ei+kei ~ a 2 DL,B-2ei+ke 2 ) /3 k (-aS L,B+ei+ke 2 + a 2 S L,B-2ei+ke 2 ) 



k£\0,l) 



2^ B k (aT UB+ei+k ? 2 - a 2 T L ^ 2ei +kt 2 ) ^ P k (- aD * »- f <! ' ' D 

te{0,l) te{0,l) 



L.B-1--+U-- ' " ^ L,B-2ei+ke 2 ) 
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Q 



BL 



X a j (fiD BtL+jei+e2 -B 2 D B x +/ei -2e 2 ) ^ a j (-/3S BMjej+ei +P S B.L+j ei -2e 2 ) 
e{0,1} j£{0,D 



M0,i} 



J£{0,D 



^ ka j p k D BtB+M+ke2 - ^ ka i p k S B ,B+j^+ke 1 

;'£{-l,0,l},te{-l,l} 7'e{-l,0,l},te{-l,l} 

,'ffl = 

j£{-l,0,l},*e{-l,l} ./e(-l,0,l|,te(-l,l) 

Finally, we discuss the B operator from (l34t . which describes the effect of the auxiliary 
densities £ on the mismatch. As with A, since the mismatch involves values on only a single 
curve <9Q, there is no opportunity for cancellation. Inserting (|4TT i into d2TT i we get 



je{0,l},te{-l,0,l} 
jE{-l,Q,l},ite{Q,l} 



DdQ,L+jei+ke 2 ~$ 
TdCl. 







L+jej+iej D *dCl,L+jei+ke 2 



DgQB+je,+ke 2 S dn,B+je l +ke 2 

TatxB +jei+kei -D* gnB+M+ke2 



(42) 



Summarizing the above, £ is a compact perturbation of the identity. Its blocks C and Q 
involve interaction distances greater than the unit cell size. Its block A involves distances con- 
trolled by the shape of the inclusion and its nearest approach to its neighboring images. Its block 
B involves distances determined by the nearest approach of 3D. to dU . 



4.3. Numerical implementation and discretization of B 

We discretize the four blocks of the integral operator E in ( [34-b to give the matrix E e 
£.(2N+4M)x(2/v+4M) as f ji ows \y e sample the densities on dQ. at equispaced points with respect to 
the given definition of the curve, as in Section |3~T1 We sample the densities on the walls L and 
B at M standard Gaussian nodes, as in Section [3721 A is then discretized in the same way as A QP 
in Section [3~Tl with a mix of the periodic trapezoidal rule and Kress' singular quadratures for the 
self-interaction of 3Q. The (Nystrom) method (l22l may be used for the off-diagonal block C, 
and also for the wall's self-interaction Q. No special singular quadratures are needed in Q, due 
to the cancellations discussed above. 

The B operator (l42l involves computing the field due to source densities on walls L and B 
(and their images shown in Fig. QJ)) at targets on dQ.. When the distance from the inclusion 
to boundary dist(3Q, dU) is large, the plain Nystrom method may be used to construct the dis- 
cretized matrix B. We will refer to this as discretization method Bl. With nodes y m and weights 
w m on wall L, and nodes Xj on <9Q, for example, the term S gn,i in the ( 1 ,2)-block of (l42l becomes 
the matrix S e C NxM with elements 8 ]m = jH^\a>\Xj - y m \)w m . 

When dist(<9Q, dU) becomes small, of course, the convergence rate of method Bl will be- 
come unacceptably poor. However, by construction, for a Bloch eigenfunction the field (PiTt 
generated by the wall densities in £ has no singularities in the 3x3 neighboring block of unit 
cells. Hence these densities remain smooth, poor convergence being merely due to inaccurate 
evaluation of their field close to the walls. This leaves room for a large number of options: 

16 



B2) For the rows of B corresponding to target points on dQ. that are distance do or closer to dll, 
use adaptive Gauss-Kronrod quadratur^l with integrand given by the product of the kernel 
function and the Lagrange polynomial interpolant 1I32I Sec. 8.1] for the density at the M 
quadrature points, do is some 0(1) constant. For the other rows, use method Bl. 

B3) Project onto an order-L cylindrical /-expansion at the origin. This is done by computing 
a representation d27b for each of the point monopole or dipole sources in the quadrature 
approximation to the source densities on the walls, and then evaluating this at the target 
quadrature points on dQ to fill the elements of B. The example term discussed for Bl gives 
S = RP, where the "source-to-local" matrix P e c( 2L+1 ) xM has elements 



P lm = l -Hf\u\y m \)e- iW '" W „ 



and converts single layer density values to /-expansion coefficients. This follows from 
Graf's addition formula yj, Eq. 9.1.79]. The expansion matrix R e C Nx ( 2L+1 ) has ele- 
ments Rji = Ji(p)\K.j\)e l1 ^' . In the above 6 m , <pj are polar angles of points y m , Xj respectively. 
Similar formulae apply for double layers and evaluation of derivatives. To reduce dy- 
namic range (hence roundoff error) we in fact scale the /-expansion by the factors p\ of 
Section l3~2l (this does not change the mathematical definition of S .) 
B4) Use a more sophisticated quadrature approach, such as those of 



Methods B2-B4 evaluate u QP due to a spectral interpolant of the discretized wall densities, with 
an accuracy that persists up to the boundary of U. Note that this does not increase the number M 
of degrees of freedom associated with each such density. Since the underlying density is smooth 
(in fact analytic), the convergence rate is high and we are able to keep M very modest. 

We have implemented methods Bl, B2 and B3. We use the quadrature weights to scale the 
matrix E to give E in an analogous fashion to (l26l i. so that singular values of E approximate 
those of E. 

Finally, there are many possible ways to locate parameter values (u>, a, b) where E is singular. 
In this paper, we will simply plot its smallest singular value cr mm (E) vs the Bloch parameters, as 
in Section[3] 

5. Results of proposed scheme 

We first test the convergence of the new scheme for the same small inclusion used in Sec- 
tion [3] with the simplest discretization method for B, namely Bl. As before, we test two Bloch 
parameter b values, one which is far from an eigenvalue, and one of which is guaranteed to be 
an eigenvalue according to Theorem |4] Fixing N = 70, which was found in Section [3~T1 to be 
fully converged when at an eigenvalue, we first vary M, the number of nodes per unit cell wall. 
Fig- shows the convergence of the minimum singular value of the discretized matrix E to its 
converged value (when far from an eigenvalue), or to zero (when at an eigenvalue). The conver- 
gence is spectral, and in both cases full machine accuracy is reached at M = 30. (For N > 70 
the results are unchanged.) Thus for a matrix of order 2N + AM = 260, we are able to locate 
the desired band structure with relative error around 10~ 15 in the Bloch parameters (a, b). Filling 



This was implemented with MATLAB's quadgk, which uses a pair of 15th and 7th order formulae, with relative 
tolerance set to 10~ 12 . 
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b) c) 




Figure 5: Convergence of new periodizing scheme using auxilliary densities (described in Section|4j, using method Bl, 
for the same geometry and parameters as in Fig. [2^. The meaning of the two curves is also the same as in the earlier 
figure, a) Absolute error in cr min (£) vs M the number of nodes on each unit cell wall, for fixed N = 10 nodes on <9fl. b) 
Same as a) except convergence vs N, for fixed M = 30. c) Same as Fig.fSp but using the new scheme: note the absence 
of pollution by the empty band structure. 



such a matrix takes around 0.45 sec and computing the complex SVD around 0.15 sec. Fur- 
thermore, by storing coefficient matrices in the expansion E = 2-i</,/t<2 a 

jpk^UM) at fixed 

oj, we 

can fill E for new a, b values in 0.05 sec. 

Fig. [5J5 shows that, with M in the new quasi-periodizing scheme sufficient to yield machine 
precision, the error convergence rate with respect to N is the same as that of the old scheme. 
Fig. [5J; demonstrates the robustness of the scheme, by plotting the smallest singular value over 
the same region of parameter space as Fig. |3j>. Notice that the location of the desired band 
structure (black line) is unchanged, but that the divergent behavior near the empty resonant band 
structure has entirely vanished. 

5.1. Inclusions approaching and intersecting the unit cell wall 

Given a crystal of inclusions, it may be impossible to choose a parallelogram unit cell U 
whose boundary does not come close to or even intersect dQ.. Although this is not an issue for 
the scheme of Section|3] for the new scheme which relies on dU it is a potential problem. 

We first show that, as expected, with method Bl the error performance deteriorates exponen- 
tially as <9Q approaches dU. In Fig. [6^ we plot the minimum singular value at a Bloch eigenvalue, 
as a function of distance d that the inclusion has been translated in the x direction (translation 
does not affect the Bloch eigenvalue.) Numerical parameters N and M are held fixed. The loga- 
rithm of the error grows roughly linearly with d and reaches 0(1) for dist(d£2, dU) - 0, indicated 
by the dotted vertical line at around d - 0.23. Method B2, also shown in Fig.|6^, uses adaptive 
quadrature for accurate evaluation of u QP in all of U. For very small d, the inclusion is still cen- 
trally located (far from the wall) and B2 is identical to Bl, with an error of 10~ 15 . The error is 
around 10" 12 as one approaches the wall (more or less independent of d), limited by the accuracy 



3 All timings are reported for a laptop running MATLAB 2008a with a 2GHz Intel Core Duo CPU. 
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Figure 6: a) Dependence of cr m[n (E) on .v-translation distance d of dSl relative to the system of Fig. [2^, for fixed N = 70, 
and M = 30. The vertical line shows where dSl starts to touch dU. B is discretized as follows: method Bl (+ symbols), 
method B2 with d = 0.2 (□ symbols), method B3 with L = 16 (o symbols), method B3 with L = 22 and M = 40 (* 
symbols). Inset shows unit cell and inclusion at d = 0.6. b) Band structure for crescent-shaped photonic crystal shown in 
c), index n = 2, shape (0.265 cos 2nt + 0.318 cos Ant, 0.53 sin 2m), < t < 1, unit cell ei = (1,0), e 2 = (0.45, 1). A tour 
rXMT of the Brillouin zone is shown, where T is (a, b) = (0,0), X is (tt, 0), and M is (n, n). d) Contours of the Bloch 
mode Re[«] with parameters shown by the dot on the band diagram. 
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of quadgk. This proves that the deterioration seen with Bl is associated with the B operator 
block, and can be remedied merely by careful discretization of B without increasing the matrix 
size. We did not bother continuing the computation with B 1 or B2 after the inclusion crosses 
the wall; here they fail because ( HTl i. as constructed, represents m qp only inside U (jump relations 
cause the values outside U to be different). We note that to use Bl or B2 correctly, one would 
have to wrap the boundary points outside U back into the cell, evaluate at the wrapped point, and 
correct for phase. Method B2 is not very useful in practice since the call to a black box adaptive 
quadrature routine causes the matrix fill time to increase to 55 sec. 

Finally we use method B3 with L = 16, M = 30 and with L = 22, M = 40. In the first 
case, errors grow slowly to around 10~ 12 as dist(<9£2, dU) reaches zero, and then continue to grow 
slowly to a plateau at around 10~ 9 , even though most of <9Q now falls outside of the unit cell. The 
cost of B3 is not much more than Bl, taking 0.7 sec to fill E. Note that the /-expansion used to 
represent m qp has effectively carried out analytic continuation beyond U. This is stable because 
our image structure has pushed the singularities out beyond the nearest image cells. It is perhaps 
worth observing that some care must be taken in setting L. With M — 30, increasing L above 
16 would worsen errors (not shown). The reason is that the coefficients \l\ > 16 involve more 
oscillatory integrands which are not resolved by M — 30 points. Increasing M to 40 permits 
increased precision with L = 22, as seen in Fig. |6^. 

There is another potential pitfall with method B3 as implemented; if both L and d get larger, 
there may arise singular values of E which become exponentially small, associated with highly 
oscillatory non-physical densities on the farthest part of dQ.. For illustration, with L = 16 and d = 
0.6, the second-smallest singular value is 10~ 4 ; with L — 22 the second smallest singular value 
shrinks to 10~ 6 . (When d = 0, the second smallest singular value is 10 _I .) This is troublesome for 
eigenvalue search methods that track cr mi „(E) vs Bloch parameters, since the desired minima will 
be obscured by these spurious small singular values everywhere except in a small neighborhood 
of the desired band structure. We will discuss search methods less sensitive to this problem in a 
future paper. For now the lesson is that, when parts of an inclusion extend far beyond U, there is 
a price to pay for making use of analytic continuation. 

5.2. Application to band structure 

We compute the band structure of a more difficult crystal in Fig. |6j5. O is far from circular, 
hence simple multipole methods i23ll would not be accurate. The closest approach to its neigh- 
bors is only 0.06, so that N = 150 points are needed in discretizing the inclusion boundary. Note 
that any parallelogram unit cell must intersect dQ, so the method of 114911 cannot be used without 
modification. We use method B3 with M — 35 and L - 18. As illustrated before in Fig.|2j5, 
the minimum values of cr Iui „(E) on the band structure indicate the size of the errors in the Bloch 
parameters found. By this measure, sampling 100 random points on the first 15 bands, we find 
a median error of 3 x 10~ 10 and a maximum 1.6 x 10~ 9 . 1.7 sec were required to fill the matrix 
E of order 440 once for a given a>, a, b (and 0.13 sec for subsequent values of a, b). The SVD 
required 0.7 sec for a matrix of this size. We located the band structure using 8000 such evalu- 
ations and a specialized search algorithm, which we will describe in a forthcoming paper. The 
search algorithm is also accelerated by computing the determinant of E rather than the SVD, at 
a cost of 0.1 sec for each matrix. The total CPU time required was 35 minutes. 

Fig-UJi shows a single Bloch mode on the 11th band for this crescent-shaped crystal. This 
took 16 sec to evaluate on a 100 x 100 grid over U using d33t . and the /-expansion for (I4T1) (with 
no fast multipole acceleration). 
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6. Conclusions 



We have presented two algorithms for locating the band structure of a two-dimensional pho- 
tonic crystal, in the z-invariant Maxwell setting. The first (Section [3]) uses the quasi-periodic 
Green's function. Theorem [4] guarantees the success of this method (no spurious or missed 
modes) as long as the band structure for the empty unit cell is avoided, where we have shown 
that the method fails. The second method (Section |4| introduces a small number of additional 
degrees of freedom on the walls to represent the periodizing part of the field: numerical evi- 
dence suggests that it is immune to breakdown for any Bloch parameters (Conjecture Q. The 
two schemes are connected by the following observation. 

Remark 8. Computing the Schur complement formula for the operator system (134-b recovers the 
quasi-periodic Green 's function approach described by ( 1201 l. In particular, 

A QP =A-BQ- l C. 

The quasi-periodic Green 's function approach fails when Q becomes singular and A QP blows up. 
The full system d34l >, on the other hand, remains well-behaved. 

We have shown spectral convergence for both schemes, achieving close to machine precision 
accuracy on simple crystals using only a few hundred degrees of freedom, hence CPU times of 
less than 1 sec for testing at a single parameter set (w, a, b). In the new scheme we have shown 
(method B3) how to handle the passage of the inclusion through the unit cell boundary, with- 
out much sacrifice in accuracy, without much extra numerical effort, and with no bookkeeping 
needed to determine which points of d£l lie in U. The latter is convenient for larger-scale or 
three-dimensional (3D) computations if existing scattering codes are to be used to fill the A op- 
erator block. Other ways to handle this intersection problem exist, such as a variant of B2 which 
wraps points on dQ. back into U, with which we have preliminary success. 

We have not discussed the methods we use for the nonlinear eigenvalue problem, due to space 
constraints. The scheme of Yuan et al [49] uses a quadratic eigenvalue problem, and factorizes 
the scattering matrix of the inclusion at each u>, hence may be faster than our scheme for small 
systems. However, moving to large-scale systems with more than 10 4 de grees of freedom, such 
a factorization would be impractical compared to an iterative version of our scheme. 

Some generalizations of what we present are straightforward, such as multiple inclusions per 
unit cell, non-simply connected inclusions, or inclusions with corners (using quadrature rules 



such as II 111 12511 ). There exist regimes, however, that would require some modification. These 
include two phase dielectrics one or more of which are connected through the bulk (sometimes 
called bicontinuous), and unit cells which are highly skew or have large aspect ratios. 

Our new representation for quasi-periodic fields can also be used for scattering calculations 
from periodic one-dimensional arrays of inclusions in 2D and one or two-dimensional arrays in 
3D. Because we rely entirely on the free-space Green's function, it should be straightforward to 
create quasi-periodic solvers from existing scattering codes. We will describe such solvers at a 
later date. 
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Appendix A. Proof of Theorem [4] 

Recall the Green's representation formulae Sec. 3.2]. If u satisfies (A + oj 2 )u — in n, 
recalling that iC and u~ signify limits on <9£2 approaching from the inside, and the normal always 
points outwards from n, then 

+ tixa < A1 » 

The exterior representation has the opposite sign: let u satisfy (A + <x> 2 )u — in K 2 \ fl and the 
Sommerfeld radiation condition, that is, 

^ - iuu = o(r l/2 ), r := |x| -> oo (A.2) 
or 

holds uniformly with respect to direction x/r. Then, 

-S^+i^V = ( ° in " - (A.3) 
y u in K \ 12 

We will need the following quasi-periodic analogues. 

Lemma 9. Let u satisfy (A + oj 2 )m = in Q., and Q. c U, Then for each Bloch phase (a,f$), 



—u in Q 
in U \ Q 



-S^u'+D^u- = { n " (A.4) 



Proof: Write G QP using dTTb and notice that each term other than (m, n) = (0, 0) contributes zero. 
This is because all points in U lie outside each closed curve <9Q - mt\ - ne2, and we may apply 
the second (extinction) case of dA.ll ) to show that they have no effect in U. □ 



Lemma 10. Let u satisfy (A + u> 2 )u = in U \ Q and quasi-periodicity ©-(O, and Q. c U. Then 

-S (W V + ^ (W V - 1 inn 
<V«„ - I u . mUxU (A.5) 



Proof: We follow the usual method of proof [ 16, Thm. 3.3] but with the quasi-periodicity con- 
dition playing the role of the radiation condition. Apply Green's 2nd identity to the functions u 
and G QP (x, ■) in the domain U \ n if x e n, or the domain {y 6 U \ n : |x - y| > s] if X 6 U \ n. 
In the latter case the limit e — > is taken, and (fTTT i shows that only the (m, ri) = (0, 0) term 
contributes to the limit of the integral over the sphere of radius e. In both cases the boundary 
integrals contain the term 



dG 

— — (x, y)M(y) - G QP (x, y)w n (y) ds y , (A.6) 
on v 



JdU 



which vanishes by cancellation on opposing walls, since u is quasi-periodic with phases (a,j8), 
but G QP (x, ■) is anti-quasiperiodic, i.e. quasi-periodic with phases (a~ l ,/3~ l ). □ 
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Turning now to Theorem [4] to prove the if part, we show that whenever the operator has a 
nontrivial nullspace, a Bloch eigenfunction u may be constructed, i.e. a solution to ©-(O that 
we must take care to show is nontrivial. Let 77 = [t; -cr] + be a nontrivial density such that 
A QP 77 = 0. Immediately we have that the resulting field u given by cfT3T > satisfies (|2]i-((9]l- We now 
define a complementary field over the whole plane minus d£l, 

= f S^cr + D^r inQ 
\ -5 (nu V - DMt inR 2 \Q 

Suppose u = 0. Then iC — u~ — and by the jump relations for S ( "^cr + T> {nui) T we get v + = — t 
andv* = cr. Similarly, since m + = u* = by the jump relations for S]£o~+ t we get v~ = — t 
and v~ = cr. It is easy to check that v solves the (swapped-wavenumber) transmission problem, 

(A + « 2 )v = inQ (A.8) 

(A + nV)v = inR 2 \Q (A.9) 

dv U7 

— intov — o(r ' ), r — > 00, uniformly in direction (A. 10) 

or 

v + -v~ = h (A.ll) 

vl-v- = W (A.12) 

with homogeneous boundary discontinuity data h - h' - 0. By uniqueness for this problem lfl6l 
Thm. 3.40] we get that v = in R 2 , from which the jump relations back to u imply cr = t = 0, 
which contradicts our assumption of nontrivial density. Thus u is a Bloch eigenfunction. 

To prove the only if part we show that, given the existence of a Bloch eigenfunction, we 
may exhibit a (nontrivial) density rj such that A QP ?7 = 0. Let w be a Bloch eigenfunction with 
eigenvalue (u>, a, b). Then let v solve ( IA.8b -( IA. 121 > with the inhomogeneous data h - -2w\dn and 
h' = — 2w n \ga.- (Note that w obeys continuity ©, hence w\gn = w + = w~ and w n \gn = = w~). 
By {uk Thm. 3.41] we know that a unique solution exists. We now claim that the densities 

o- = w n \m + v+ (A. 13) 

r = -w\ dn -v + (A. 14) 

generate precisely the eigenfunction w, i.e. the representation u of ( fl"5l > obeys u = w in U. 
We show this by substituting the densities into $15[ , then applying iA.li and ( IA.3b in £2, and 
Lemma[10]in 17 \ Q: 

f ^""V.bn - £> (n0J) w\dn + 5 (nw) v+ - £) M v + inO 

w in Q 

-w + S^v* - D^V mU\H 

On the remaining term, we use v's known jumps h and h' to get 

<-V v « _ -^qp v - *V v n ~ -°qp v - ^qp Wnlan + 2iJ Q p wlan 
= — 2w 
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where we applied Lemma[9]to the first pair, and Lemma[TO]to the second as before. Substituting 
this above shows that u = w in U. Since w has zero mismatch, the density vector rj := [r; -<t] 
satisfies A Q? j] = 0. Finally, rj must be nontrivial since rj = would imply u = by ( fT3T > which 
contradicts it being equal to the eigenfunction w. □ 

We close with a couple of remarks about the proof. Barring a sign, v in ( IA.71 > is the extension 
of u's representation ( TT3T > into its nonphysical regions, a trick originating, in the homogeneous 
context, with the proof in [16, Thm. 3.41]. Because ( fT31 l uses G QP outside, but G inside, the 
complementary problem is a nonperiodic transmission problem, which has known existence and 
uniqueness. The related analysis of i45ll uses G QP both inside and outside. This results in a 
periodic problem as the complementary problem, and it is not so clear that one can eliminate the 
possibility of spurious modes. 
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